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ABSTRACT 

We study the covariance matrix of the cluster mass function in cosmology. We adopt a two-line 
attack: firstly, we employ the counts-in-cells framework to derive an analytic expression for the co- 
variance of the mass function. Secondly, we use a large ensemble of iV-body simulations in the ACDM 
framework to test this. Our theoretical results show that the covariance can be written as the sum 
of two terms: a Poisson term, which dominates in the limit of rare clusters; and a sample vari- 
ance term, which dom inates for more abundant clusters. Our expressions are analogous to those of 
IHu fe Kravtsovl (|2003l ) for multiple cells and a single mass tracer. Calculating the covariance depends 
on: the mass function and bias of clusters, and the variance of mass fluctuations within the survey 
volume. The predictions show that there is a strong bin-to-bin covariance between measurements. 
In terms of the cross-correlation coefficient, we find r > 0.5 for haloes with M < 3 x lQ 14 h- x M Q at 
z = 0. Comparison of these predictions with estimates from simulations shows excellent agreement. 
We use the Fisher matrix formalism to explore the cosmological information content of the counts. 
We co mpare the Poisson likelihood model, with the more realistic likelihood model of iLima fc Hul 
(2004), and all terms entering the Fisher matrices are evaluated using the simulations. We find that 
the Poisson approximation should only be used for the rarest objects, M > 3 x 1O 14 /i _1 M0, otherwise 
the information content of a survey of size V ~ 13.5 h~ 3 Gpc 3 would be overestimated, resulting in 
errors that are ^2 times smaller. As an auxiliary result, we show that the bias of clusters, obtained 
from the cluster-mass cross- variance, is linear on scales > 50 ft Mpc, whereas that obtained from the 
auto- variance is nonlinear. 



1. INTRODUCTION 

The last decade of research in cosmology has largely 
been focused on devising probes to reveal the physical 
nature of dark energy and the origin of the accelerated 
expansion of the Universe. Among the most promis- 
ing pr obes, as identified for example in lAlbrecht et al.l 
(200(3), are cluster counts. 

From a theoretical perspective, the abundance of clus- 
ters per unit solid angle rffi, is an integral of the mass 
function over mass M and volume element dV: 



dN 

~dn 



dz- 



dV 
dfldz 



dMn(M) , 



(1) 



M th (z) 



where M t h{z) is a redshift-dependent mass detection 
threshold for the clusters and where the mass function is 
defined as the number of halos per unit volume and unit 
mass, i.e. n(M) = dN/dV/dM, with M the virial mass. 
For a wide range of cosmological models, n(M) can be 
accurately predicted from the semi-analytical prescrip- 
tion s based on the spherical or ellipsoidal collapse m odel 
(e.g. lPress fe Schechter|[T974l : ISheth fe Tormenlfl999l and 
see £|3.4I for more details.). 

The mass function is primarily sensitive to the statis- 
tics of the initial conditions and to the amplitude as and 
shape of the matter power spectrum; which in turn de- 
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pends on the matter density of the Universe Q m , the 
Hubble parameter h, the spectral index of the primor- 
dial power spectrum n; and the dark energy equation of 
state w = Pw/Pw, where P w and p w are the pressure and 
energy density of the dark energy. The volume element 
integral in the above equation renders the cluster counts 
even more sensitive to VL m and w. Measuring the cluster 
abundance at different redshifts can constrain a dynam- 
ical w and thus enable one to differentiate between a 
cosmological constant A and alternative dark energy sce - 
narios such as quintessence (jWang fe Steinhardtl 11998). 
or dark energy inhomq geneities coupling to dark matter 
(M anera fc Motall2006jL 

For many decades the study of clusters of galaxies has 
been a centerpiece for observational cosmology, which has 
produced many important results and cosmological infer- 
ences. Currently there are four observati onal strategies 
for de tecting clusters: A-ray emis sion (seelBorgani et all 
20011 iReiprich fc Bohringeril2002t [Schuecker et al.l 120031: 



Allen et al.l l20qf IHenrvl l200i 
Vik hlinin et al.l 120091: iMantz et aT 



Mantz et al.l 120081: 
. 120101 and refer- 
ences therein); optical emission (see Idadders et aD 
l2007t iRozo et all I2010L a nd references therein): the 
Sunyaev-Zel'Dovich effect (jSunvaev fc Zeldovichl I1972L 
hereafter SZ effect), i.e. the up-scattering of 
CMB photons off hot electrons in the intraclus- 
ter medium (see iVanderlinde fc The SPT Collaboration 



2010; Planck Collabora tion et al.1 120111: Muchovej et al 
20L1 ISehgal fc The ACT Collaborationl 120 111 and ref- 



erences therein); weak gravitational lensing (see 
iSchirmer et "al|[2007t lAbate et al.|[200l llsrael et aJ J20lol 
and references therein). 

One of the most challenging aspects of deriving cos- 
mological constraints from cluster counts is the fact that 
virial masses are not directly observable: a conversion 
is needed to translate observables such as flux, lumi- 
nosity, temperature and SZ decrement into mass. The 
mass-observable relation is de generate with cosm olog- 
ical parameters, as shown in iLima fe Hul ([2005), and 
can severely degrade the inferred constraints. Substan- 
tial progress has been made in calibrating the mass- 
observable relation in the recent years, through nu- 
merical simulation s, or by comparing diffe r ent methods 
again s t each other (iZhang et al.ll2007U2008t lOkabe et all 
l2010f ). ILima fc Hul (|2005h also proposed a self-calibration 
technique that uses the clustering of clusters to break 
the degeneracy between the uncertainties in the mass- 
observable relation and cosmological parameters. 

Owing to observational challenges, the cluster studies 
mentioned earlier employ small numbers of massive clus- 
ters (at most a few hundreds, but in general a few tens) 
to constrain cosmology. In obtaining these constraints 
it is widely assumed that the likelihood function for the 
selected clusters follows the Poisson distribution. Whilst 
this assumption may be reasonable for the most massive 
clusters, M ~ 10 15 ft. _1 M©, it will certainl y fail at lower 
mass e s. Future surv eys, suc h as eROSITA (iPredehl et all 
MHh , ESS 3 jSOp, E uclid (IRefregier et al.ll2010D . Pa~ 
STARRS3,|DH (HCKl), will be able to detect large sam- 
ples of intermediate-mass clusters, M ~ 10 14 /i _1 M©. In 
order to make accurate inferences from this data, the 
cluster likelihood function will require a more complex 
statistical treatment, and in particular knowledge about 
the covariance matrix of the mass function. 

This paper is driven by the following two questions: 
What is the covariance matrix for measurements of the 
mass function? How much are forecasted errors, which 
rely on the Poisson approximation, affected by more re- 
alistic modelling of the cluster likelihood function? The 
main theoretical tools that we shall employ to answer 
these questio ns will be the counts-in-cell formalism in- 
troduced by iPeeblesl (fl980) and further d eveloped by 
IHu fc Kravtsovl (|2003l hereafter HK03) and ILima fc Hul 
( 2004 hereafter LH04). We shall also compare the the- 
oretical results obtained via this formalism to measure- 
ments obtained from a large ensemble of A^-body simu- 
lations. 

As this pap e r was nearing submission, a study by 
IValaeeas et aLl (|201lD was reported. This work explores 
related, but complimentary, questions to those presented 
here. 

The paper is structured in the following way: in ij2]we 
review the counts-in-cells formalism and also the exten- 
sion to the cluster likelihood developed by LH04; in S}3]we 
derive the mass function covariance in a formal way; in 
21 we describe the numerical simulations from which we 
measure the mass function covariance; in $5] we present a 
comparison between the measured and the predicted co- 
variance, and in Sj6] we use the Fisher-matrix formalism 
to estimate the impact that the full covariance matrix of 

1 http:/ /pan-starrs. ifa. hawaii.edu 



the mass function has on cosmological constraints. Fi- 
nally, in <J7]we discuss and summarize our findings. 

2. THEORETICAL BACKGROUND 
2.1. The cellular model 

In this section, we give a short description of the 
counts-in-cell formalism, used by HK03 to compute the 
linear-theory sample variance of cluster counts, and by 
LH04 to estimate the impact of the latter on Fisher ma- 
trix predictions. 

Consider some large cubical patch of the Universe, of 
volume V^, and containing A~ clusters that possess some 
distribution of masses. Let us subdivide this volume into 
a set of N c equal cubical cells and the mass distribution 
into a set of N m mass bins. Let the number of clusters 
in the i th cell and in the a th mass bin be denoted Aj jQ! . 
We shall assume that the probability that the i th cell 
contains Ni ia clusters in the mass bin a, is a Poisson 
process: 



Ni 



P(N i;a \ 



exp(-m iiQ ) 



(2) 



For any quantity X, we denote the average over the 
sampling distribution-the Poisson process in this case-as 
(X) p , and the ensemble average over many realizations 
of the density field as X = (X) , termed sample variance 
in HK03. The average of Ni a over the sampling distri- 
bution can be wri tten as (see also iCole fe Kaiser! 119891 : 
IMo fc Whit3H99a >: 



, Q [l + b a 5v(xi)] , 



(3) 



where 



n a Vi is the ensemble- and Poisson- 



averaged number of counts in cell i and mass bin a. The 
volume of the cell and the cell-averaged overdensity are 
given by, 

Vi = J d 3 xiy(x| Xi ) ; (4) 

S v {x i ) = yJd a xW(x\x i )S{x) . (5) 

where W(x|x$) is the window function for the ith cell (see 
£13.21 for more details). The number density and linear 
bias of the clusters averaged over the mass bin a are 
given by: 



M a +AM a /2 

M a -AM a /2 
2 r M a +AM a /2 



b a = — 

n a JM a -AM a /2 



dMn(M) ; (6) 
dMb{M)n(M) , (7) 



where b(M) is the linear bias of haloes of mass M. 

As was shown in HK03, the correlations in the under- 
lying density field induce a correlation in the number 
counts of the cells, defined as: 



S' 



a/3 . 



: ((Ni a - m i>a ) {Nj t p - m j<f3 )) p s 

/d 3 k 
—^*(k)^-(k)P(fe) , (8) 



where for independent Poisson processes the probability 
P(N ita ,Nj t p\mi !a ,mj t p) = P{N ita \m ita )P{Nj i p\mj t p). 
In the last line we introduced the power spectrum P(k) 
as the Fourier transform of the correlation function £, 



d 3 k 
(2^)3 



P(fc) exp (-• ik • r) . (9) 



Wi(k) is the Fourier transform of the cell window func- 



tion and r 



(see 



2.2. The Gauss-Poisson likelihood function for counts 

in cells 

The likelihood of drawing a particular set of clus- 
ter counts N e {N\,\, . . . ,Nn c ,i, Ni,2, ■ ■ ■ ,N Nc , Nm } in 
the cells, given a model for the counts in the cells 
Tn g {mi i , . . . , rnN c i, rnjv c % . . . , mjv c jv m } was written 
byLH04: 



£(N|m,S) 



i.a TTli.a J 



.a— 1 i — l 



G(m|m,S) 
(10) 



with AT = N c x iV m , and where it was assumed that the 
statistics of the cell-averaged density field are described 
by a multivariate Gaussian: 



G(m|m,S 



-N/2 



IS] 1 / 2 



■ exp 



2 V 



m) T S" 1 (m 



with S defined in Eq. 



(11) 

Measurements of the bis- 
pectrum of the CMB have shown that the statistics 
of the initial fluctua tions are very nearly Gaussian 
(jKomatsu et all 12010( 1 . Whilst we know that nonlinear 
growth of structure in the present epoch drives the statis- 
tics of the density field to become non-Gaussian, in the 
limit that the cells are large compared to the coherence 
length of the field, we expect that the Gaussian approx- 
imation will be very good. 

At this point we may also be more precise about what 
we mean by ensemble and Poisson averages: 



(X(N)) p<s = J2 •■• E £(N|m,S)X(N). 

(12) 



ATi ( i=0 N Nc ,N m =0 



Equation (| 10[) can be simplified in two limits: 

• Case I: In the limit that the ensemble average vari- 
ance is much smaller than the Poisson variance: 
i.e. Su <C rn~i. In this case, the Gaussian effec- 
tively becomes a delta function centred on m and 
the likelihood simply becomes a product of Poisson 
probabilities: 



N m N c 



£(N|m)« JJ Y[P(N^ a \m^ a ) 



(13) 



a—l i—l 



• Case II: In the limit that the number of counts in 
each cell and mass bin is large, then the Poisson 
process becomes a Gaussian: 



N m N c 



Et Yl p (Ni, a \mi, a ) « G(N|m, M) 



(14) 



where M —> = 8^6% pm^a.- Hence, as shown 
in LH04, the likelihood function becomes, 



£(N|m,S) 



d A/ mG(N|m,M)G(m|m,S) (15) 



and via the convolution theorem this can be ap- 
proximated as a Gaussian with shifted mean and 
augmented covariance matrix: 

£(N|m,S) w G(N|m,C) ; C = M + S, (16) 

where M -> = s *j S */3 m i,of Note that in the 
above equation, the approximate sign is used since 
negative number counts are formally fo rbidden (for 
a mo re detailed discussion of this see IHu & Cohnl 



3. COVARIANCE OF THE MASS FUNCTION 

The final result of SjH is that in the limit of a large 
number of counts per cell, the joint likelihood for all the 
cells is a Gaussian with model mean m and with a co- 
variance matrix, C = M + S. In the following section, 
we shall use these results to answer the question: What 
is the covariance matrix for measurements of the mass 
function? 

3.1. A formal approach 

The mass function n(M) is the number density of clus- 
ters in a volume V , per unit mass. Using our counts in 
cells distribution, an estimator for the mass function in 
the i th cell is, 



hi(M a ) 



ViAMa ' 



(17) 



which, if we average over all cells and all cells have equal 
volume, becomes 



(18) 



The above estimate is unbiased, and its expectation value 
n(M a ) = (n(M a )) Ps can be formally calculated using 

Eq.(2): 

1 v i, a 



n(M a )= £ ... Yl £(N|m,S)]T 

Ni .i=0 NN c ,N m =0 i 



VuAM a 



/oo 
d^mG(m|m,S) ^ P(iVi,i|mi,i) . . . 



N, i=0 



N 



E P(NN c ,N m \m NctNm )J2 



Ni. 



N\ 



V^AM a 



V, 

N c 

E m i,a 

2 — 1 ^ 



(19) 



In a similar fashion, the covariance matrix of the cluster 
mass function can also be calculated: 



M a0 = ([n(M a ) - n{M a )\ [n{M p ) - n(M^)]) 



a—l i—l 



s,P 



4 



___ (N ia N i0 ) p 
= Y V 9 3 ' Ns ' P ~n(M a )n(Ms), (20) 
^V 2 AM a AM p { aJ K Ph V ' 

where the expectation of the product of the counts can 
be written 

oo oo 

E(W,rI] ■•• E AN TTT.Si^.V,..,.V,. 

*,3 iVi,i=0 JVjv e ,N m =0 ij" 



y d A/ "™G(m|m, S) 



E fm,«m jt p + ( N i, a ) 



(21) 



Recall that mi, a = n(M a )AM a Vi and that for the Pois- 
son distribution we have: {X 2 ) = (X) [1 + (X)]. On in- 
serting these relations into the above equation, and on 
completing the sums, we find: 

( N i,<xN jlf i} s ^[ d"mG(m\m, S) 

X 2_s [ m h<x m j,0 + m i,a § £j § a,ff] 

= E [ S if + m i,cMj,ff +' ff H,aS^jS„ t p $2) 
ij 

where in the last line we used Eq. ([5} . On inserting this 
result back into Eq. (|20[) , we obtain 



M a p = Y 



fif 777 



y/AM a AM^ 



= Qg4) n(M a )n(Mp)b a bp 
V^AM a V 2 

/ 70W;(k)^(k)F(fc) . (23) 

Considering the first term in the above, we may simplify 
this expression by performing the sums over i and j, and 
the window functions i.e. 

Y Wi(k) = 2 V5 / d 3 xexp [ik • x] VF(x|x 4 ) 

i i ^ 

= /" d 3 xcxp [ik • x] ^ W(x|xj) = V„W{k). 



(24) 



Hence, we have that the covariance matrix can be writ- 
ten: 



M afj = n{M a )n{Mp)b a bpa*(V») + 



VnAM a 



-, (25) 



where <7 2 (V^) is the mass density variance in the entire 
volume 



0K(k)| 2 P( fc) 



(26) 



From Eq. (|25|) it can be seen that the crucial quan- 
tity which controls the covariance between estimates of 



the mass function in different mass bins is a(V^). The 
strength of the covariance is also modulated by the linear 
bias and the mass function in each of the bins considered. 

3.2. A short-cut to the covariance 

Whilst in the above we have presented a formal deriva- 
tion of the mass function covariance from the HK03 and 
LH04 formalism, there is a more intuitive approach to 
arriving at the same result as given by Eq. (I25|) . which 
we now mention. 

Let us consider the limiting case where we have a single 
cell that fills the whole of our sample space Vi — > V^; 
also rtii^a — ¥ m a and similar for all the other quantities 
defined in the cells. The above formalism still applies, 
and we have that the covariance matrix of mass function 
can be written: 



M a f3- 



S 



V 2 AM a AM« 



■-if "H,a 

a ^V 2 AM a AM« 



-ri{M a )ri(Mf i )b a b si a' 2 {V ll ) + 8, 



n(M a ) 



where <r(V| J ) is the variance in the total volume. 

3.3. The cross-correlation coefficient 

As a direct corollary to the previous results, we may 
write an expression for the correlation matrix, which is 
defined 

r a s = ■ . (28) 

On factoring out [nfMj/^AMj 1 / 2 from y/M aa in the 

denominator, and a similar term from v M.W , and on 
using the fact that m a = n(M a )AM a V^, we find 



r aff 



y/m a mf3 babpa 2 ^) 



-2 



1/2 



1 + m a b a a 2 (V ll ) 1 + m / 36 /3 cr 2 (V M ) 



1/2 



____ ( (29) 
Two limits are apparent: when ^m a mp babpa 2 ^^) <C 
1, then r a fj — > 5^ a and the mass function covariance 
matrix is decorrelated; this would happen for the case of 
rare halos, for which the mass function is very small. On 
the other hand, when y/m a Wip b a bi3<j 2 (V^) S> 1, then 
faff 1 and the covariance matrix is fully correlated. 
This would be the case for smaller halos, for which the 
mass function is quite large. 



Finally, we note that taking 



and hence 



o~(Vu) — >■ 0, does not guarantee that the correlation be- 
tween different mass bins is negligible. As the above 
clearly shows, it is the quantity V^<r 2 (V^) that is required 
to vanish for negligible correlation to occur. For a power- 
law power spectrum, we would have that V^a 2 (Vff) oc 
i? 3 i?~( 3+ ") oc R~ n , which only vanishes for n > 0. For 
CDM we have a rolling spectral index, and n > for k < 
0.01 /iMpc~\ which implies that L box > 500/i _1 Mpc for 
the covariance to diminish. 

At this juncture, we point out that Eqs ([23)1 and (|2"9"j) 
constitute the main analytic results of this work, and all 
which follows will be concerned with their validation and 
implications. 
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TABLE 1 

zHORIZON COSMOLOGICAL PARAMETERS. COLUMNS ARE: DENSITY PARAMETERS FOR MATTER, DARK ENERGY AND BARYONS; THE EQUATION 
OF STATE PARAMETER FOR THE DARK ENERGY W; NORMALIZATION AND PRIMORDIAL SPECTRAL INDEX OF THE POWER SPECTRUM; 

DIMENSIONLESS HUBBLE PARAMETER. 



Cosmological parameters 


S~2m 


Qde 


n b 


to 


o"8 


n 


J?o[kms- 1 Mpc- i ] 


zHORIZON-I 


0.25 


0.75 


0.04 


-1 


0.8 


1.0 


70.0 


zHORIZON-Vla/Vlb 


0.25 


0.75 


0.04 


-1 


0.8 


0.95/1.05 


70.0 


zH0RIZ0N-V2a/V2b 


0.25 


0.75 


0.04 


-1 


0.7/0.9 


1.0 


70.0 


zH0RIZ0N-V3a/V3b 


0.2/0.3 


0.7 


0.04 


-1 


0.8 


1.0 


70.0 


zH0RIZ0N-V4a/V4b 


0.25 


0.8 


0.04 


-1.2/-0.8 


0.8 


1.0 


70.0 



TABLE 2 

zHORIZON NUMERICAL PARAMETERS. COLUMNS ARE: NUMBER OF PARTICLES, BOX SIZE, PARTICLE MASS, FORCE SOFTENING, 

NUMBER OF REALIZATIONS, AND TOTAL SIMULATED VOLUME. 



Simulation Parameters 


-/Vpart 


L sim [Mpch- 1 ] 


m p [/i- 1 M ] 


l soit [kpc/i -1 ] 


■^V ensemta 


Vtotl^Gpc 3 ] 


zHORIZON-I 


750 3 


1500 


5.55 x 10 11 


60 


40 


135 


zHORIZON-Vl, -V2, -V4 


750 3 


1500 


5.55 x 10 11 


60 


4 


13.5 


zHDRIZQN-V3a 


750 3 


1500 


4.44 x 10 11 


60 


4 


13.5 


zHDRIZDN-V3b 


750 3 


1500 


6.66 x 10 11 


60 


4 


13.5 



3.4. Ingredients for evaluating the covariance 

To evaluate the covariance matrix we need to provide 
models for n(M), b(M) and the Fourier transform of the 
survey window function. 

To compute n(M) and b(M) we employ the mass 
functi on and bias models presented in iSheth fc Tormenl 
{USED: 

dn p dlogf 

dl^M ~ M /sT( ^dbiM 5 (30) 



f ST (v)=AJ^-v[l 



6st(^) = 1 



qv 



{qv 2 ) p ] exp 

2p/S, 



qy 

2 



- 1 



Shc. 



1 + {qv 2 )P ' 



(31) 



(32) 



where A = 0.3222, q = 0.707, p = 0.3. In the above 
we have introduced the peak-height v(M) = S sc /a(M), 
where 5 SC — 1. 6&6/D(z) is the spherical overdensity for 
collapse, and where er 2 (M) is the variance of the linear 
density field extrapolated to z = 0, smoothed with a 
spherical top-hat filter of radius R (see below for more 
details). This radius is defined so as to enclose a mass 
M = 47rpi? 3 /3, with p the mean matter density of the 
Universe at the present epoch. 

For the survey window function we shall consider two 
simple examples. The first is a cubical top-hat, defined 
by: 



0, 



ibox/2 < X 1 < x\ + ibox/2 



otherwise 



(33) 

where / € {1,2,3} denotes the Cartesian components of 
the vectors, j is the cell index, and Lt,ox i s ^ ne s i ze °f ^ ne 
cell of volume V, = L^ m . The Fourier transform of this 



top-hat window function is: 



Wj(k) — exp(ik • x 3 



1=1 



(34) 



where jo(y) = smy/y is the zeroth order spherical Bessel 
function. The volume variance for this window function 



is: 



(A;i,fc 2 ,fc 3 )|iy(k)| 2 , 



sn <; / ^^(fc 1 ,fc 2 ,fc 3 )iM/(k)i 2 (35) 



1=1 



2tt 



where in the second equality we have used the isotropy of 
the power spectrum, e.g. P(kx,k2,ks) = P(—ki,k2, k^). 
In Eq. (1331) we use the following relation: 



r(k)i 2 =n 



1=1 



kiL hc 



(36) 



The second window function is a spherical top-hat: 
/3/(4^ J R 3 ),|x J | <r< \ Xj \ + R 



Wj(r) 



\0, 



otherwise , 



where R is the radius of the spherical top-hat. The vari- 
ance of the density field in this case has the familiar form: 



1 



dkk 2 P(k)W 2 (kR). 



for which the Fourier transform is given by 
3 



W(x) 



SIM — XCOSX] 



kR . 



(37) 



(38) 



On a technical note, we point out that for the k- 
space integrals given by Eqs (l3"5j) and (|3"T|) , we have in- 
troduced lower and upper limits k m - m > and fc max , re- 
spectively. For a real survey, the upper limit is decided 
by the resolution of the instrument used. If the mea- 
surements are made from numerical simulations, which 
is the case with this work, the softening length of the 
simulations will dictate the largest frequency Fourier 
mode available: fc max = 27r// so f t and for our simulations 
fcmax ~ 100 h Mpc -1 . However, in practice the largest 
useful Fourier mode is much smaller, and occurs where 
the shot-noise correction to t he power spectrum becomes 
comparable with the signal (|Smith et alj r2003). 
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The lower limit k m - m is a more complex issue. In the 
case of simulations, no modes with wavelength larger 
than the simulation box L s i m can contribute to the 
variance, which suggests the straightforward solution of 
adopting fc m ; n = 2n/L s i m . Since we are attempting to 
confront the theory with the reality defined by simula- 
tions, we shall always assume this cut-off scale. How- 
ever, for real surveys, the variance on a given scale will 
be affected by the existence of modes on scales larger 
than the size of the survey. We therefore recommend in 
this case fc m i n — > 0, or at least the inverse horizon size 
at the redshift of the survey. For more discussion of the 
importance of fc m i n for the predictions of the variance, 
see discussion in Appendix [AJ 

Note that in the above we shall relate the radius R of 
the spherical top-hat to that of the cubical top-hat func- 

1/3 

tion, through the relation R — (3/47r) ' Lbox- In other 
words the volumes of the spherical and cubical sample 
volumes are taken to be identical. 

4. TV-BODY SIMULATIONS 

We study the covariance matrix with a suite of 40 
large numerical simulations, executed on the zBOX-2 and 
zBOX-3 supercomputers at the Institute for Theoretical 
Physics, University of Zurich. For all realizations snap- 
shots were output at: z = {5,4,3,2,1,0.5,0}. We shall 
refer to these simulations as the zHORIZON Simulations. 

Each of the zHORIZON simulations wa s performed us - 
ing the publicly available Gadget-2 code (Springel 2005), 
and followed the nonlinear evolution under gravity of 
N = 750 3 equal-mass particles in a comoving cube of 
length L s ; m = 1500 /i _1 Mpc. The cosmological model 
is sim ilar to that determin ed by the WMAP experi- 
ment (|Komatsu et al.l [2009) . We refer to this cosmol- 
ogy as the fiducial model. The transfer function for the 
simulations was generated using the publ i cly available 
cmbf ast code (jSeliak &; Zal darriaga 1996; Selia k~et~aT] 
2003), with high sampling of the spatial frequencies on 
large scales. Initial conditions were set at redshift z = 50 
using the serial version of the publicly ava ilable 2LPT 
code (|Scoccimarrolll998t ICrocce et ail 120061) . Table □ 
summarizes the cosmological parameters that we sim- 
ulate and Table [5] summarizes the numerical parameters 
used. 

In this paper we also study the Fisher matrix of 
cluster counts for which we use another series of sim- 
ulations. Each of the new set is identical in every 
way to the fiducial model, except that we have varied 
one of the cosmological parameters by a small amount. 
For each new set we have generated 4 simulations, 
matching the random realization of the initial Gaus- 
sian field with the corresponding one from the fiducial 
model. The four parameter variations that we con- 
sider are {n {0.95, 1.05}, a 8 {0.7,0.9}, n m -> 
{0.2,0.3}, w -> {-1.2,-0.8}}, and we refer to each 
of the sets as zH0RIZ0N-Vla,b,. . . ,zH0RIZ0N-V4a,b, re- 
spectively. Again, the full details are summarized in Ta- 
blesffl&H 

Lastly, dark matter halo catalogues were generated 
for all snapshots of each si mulation using th e Friends- 
of-Friends (FoF) algorithm (Da vis et al.|[l985 f). with the 
standard linking-length parameter b = 0.2, where b is 
the fraction of the inter-particle spacing. For this we 
employed the fast parallel B-FoF code, kindly provided 



to us by V. Springel. The minimum number of particles 
for which an object is considered to be a bound halo was 
set at 20 particles. This gave a minimum host halo mass 
of M ~ 1O 13 M //i. 

5. RESULTS 

In this section we confront the counts-in-cells theory 
with the results from TV-body simulations. 

5.1. Cell variance in simulations and theory 

Since tr 2 (V^) plays a vital role in determining the 
strength of any covariance in the mass function measure- 
ments, we shall make a detailed study of it, for both win- 
dow functions discussed in Sj3]and considering volumes of 
varying size. We evaluate tr 2 (V r At ) in two different ways, 
analytically and from iV-body simulations. Furthermore, 
owing to concerns regarding the impact of nonlinear bias 
and mass evolution, we also compute the matter-matter, 
halo-matter, and halo-halo variance, which we denote as 
fmm(^), a hm( V ») and CT hh(^)> respectively. Compar- 
ing these quantities will then make clear any departures 
from linearity. 

Our analytical approach to determining the variances 
is based on standard quadrature routines to evaluate the 
theoretical expressions: for Eq. (|3"5)) . we use the multi- 
dimensional Monte-Carlo integration routine VEGAS; and 
for Eq. (l37l). we use t he QR0MB routine (for more details 
see iPress et all Il992| ) . In evaluating these integrals we 
take the linear theory power spectrum matching our sim- 
ulations, fully described in $4l Also, we take the largest 
mode in the simulation box to determine the lower limit 
of the fc-integrals fc m ; n . 

The second method is one of brute force: we measure 
fmm^): ^hmt^) and °hh(%) directly from the ensem- 
ble of simulations. Our estimator for the variances can 
be expressed as: 

f d 3 k 

4> = J J^Pab(k)W 2 (kL hox ) 

N s /2 

w^r J2 Ab(kiifc)|W(fc«fc,Lbo«)| 2 ,(39) 

M i,j,k=-N g /2+l 

where the indices k) label the Fourier mesh cell and 
kijk the magnitude of the wavenumber corresponding 
to that cell. The total number of grid cells considered 
is iV|; also, a and b are <E {m, h}, and P a &(kyfc) = 
V r A1 <5a(kyfc)*5b(kyA : ) are estimates of the various auto- 
and cross-power spectra. The window functions are as 
given in The estimates of the variance also require a 
correction for shot-noise, which for the halo-halo variance 
we implement in the following way: 

<^h,c = ^hh,d - jjr E \W(ki jk ,L box )\ 2 , (40) 

where Nh is the number of halos in the considered mass 
bin, and o~wu c and <r 2 h d are the variance of the continu- 
ous and discrete halo density fields, respectively. There 
is a similar shot-noise correction for the matter-matter 
variance; we assume that the halo-mass cross-variance 
requires no such correction. Note that the above method 
for estimating <t(VJj) is not the conventional one, where 
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Fig. 1. — The r.m.s. density variance as a function of the 
sample volume size Lbox- From top to bottom, we show re- 
sults for <r mm (V^), a hm (V fl ), o-hh(V M ), respectively. In each 
panel, blue empty and solid red circles denote measurements 
from the simulations, made using the spherical and cubical 
top-hat filter functions. The corresponding analytical predic- 
tions for the variance are denoted by the dashed blue and 
solid red lines, respectively. The size of the simulation box 
1500/i _1 Mpc is indicated by a black vertical lines, and the 
measurements are an average of 40 simulations. 



one partitions the real space counts into cells and then 
computes the variance of that distribution. However, it 
should be entirely equivalent, but with the added ad- 
vantages of being fast, since we are using an FFT, and 
allowing for the computation of the variance in arbitrary 
cell structures. 

Rather than testing all of the halo mass bins that we 
will employ later for the mass function covariance, we 
have chosen to show results for all the haloes in the sim- 
ulation with M > 1O 13 /i _1 M0. Figure Q] presents our 
results for CT mm (V^), a hm (V r M ) and (Thh(V^) as a func- 
tion of the cubical window function size, Lbox! recall 
that for the spherical window we take the radius to be 

1/3 

B = (3/47r) ' Lbox- In all three panels, the points rep- 
resent results from the iV-body simulations, whereas the 
lines denote the analytical integrals. The red full circles 
and solid lines are obtained by smoothing the density 
field with the cubical top-hat, while the blue empty cir- 
cles and dashed lines denote smoothing with the spher- 
ical top-hat function. The simulation results represent 
the mean of the 40 realizations, with errors appropri- 
ate for a single run. The size of the simulation box 
(L s ; m = 1500 /i~ 1 Mpc) is marked through a vertical black 
line on the horizontal axis. The effects of the shot-noise 
corrections on the estimates of d^ amc and <r^ h c are too 
small to be noticed on this log-log plot. 

As expected for a hierarchical mass distribution, in all 
cases the variance decreases steeply with the increasing 
box size. On comparing the results obtained from the 
simulations for the two window functions, we find very 
good agreement up until the size of the cubical region 
becomes similar to the size of the simulation cube. At 
this scale, the variance from the cubical window func- 
tion displays a significant loss in signal. For scales larger 
than the simulation box, the smoothing result become 
somewhat meaningless and unstable due to the oscilla- 
tory nature of both window functions, which can be seen 
from the measurements. 

Turning to the evaluation of the theoretical expressions 
for the variance, we see that, in the case of the spherical 
top-hat there is excellent agreement between the simula- 
tions and the theory on small scales, Lbox < 200 h~ 1 Mpc. 
For Lbox > 200/i _1 Mpc, the linear expressions underes- 
timate the measurements by rs 20% or even more. How- 
ever, on comparing the theoretical predictions for the 
cubical filter function with the measurements, we find a 
large discrepancy. We tested whether this was due to 
an error in the VEGAS evaluation of the integrals. An in- 
dependent check with mathematica produced the same 
results. 

After some investigations, we found that the discrep- 
ancy between the simulation and theory results was 
solely attributable to the difference between the discrete 
lattice structure of the Fourier space used in the simula- 
tions, and the continuum of Fourier modes used in the 
numerical integrals. A detailed discussion of this is pre- 
sented in Appendix IA.ll In that section we also show 
that as the simulation box size is increased, the theory 
and simulation results converge. Further, as is shown in 
Appendix I A . 2 1 the theory predictions are sensitive to the 
lower limit fc m i n . In applying this to the real Universe, 
we suggest letting fc min — > 0. 
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Fig. 2. — Comparison between the halo bias measured from 
the simulations and the Sheth-Tormen linear theory predic- 
tions as a function of the sample volume length. The symbols 
are as in the previous figure, and the theoretical prediction is 
represented by the dashed green line. The top panel shows 
the bias derived from the halo-matter variance, while the bot- 
tom panel shows the bias from the halo-halo variance. The 
lower panel also shows the importance of the shot-noise cor- 
rection on the &hh measurements: the upper and lower sets 
of points denote the halo-halo bias before and after the shot- 
noise correction, respectively. 



5.2. Linearity of the bias 

In linear theory, the relation between the variances 
plotted in Figure [1] is given by: 

= boUyj = 1 ^jy„) . (4i) 

b is the average linear bias from Eq. (JT)), estimated 
for a single mass bin containing all halos larger than 
10 13 h~ 1 MQ. F or the theoretical bias, we use the Sheth- 
Tormen model ISheth fc Tormenl (fl999l . |2003) . averaged 
over the same mass bin. All quantities are at redshift 
0. Since the bias is > 1, Chh(V AI ) is slightly larger 
than UhmCV^), which in turn is also slightly larger than 
cr mm (V^). At this level of detail the differences between 
the curves appear to be well related to each other as in 
Eq.flU]). 

To check this more accurately we next estimate the 
halo bias in the simulations and compare it directly with 
the theoretical predictions. In direct analogy with the 
Fourier-space bias estimates in iSmith et al.l (|2007l) . we 



construct the following real-space bias estimates: 



b hm = 4 s1 ; &hh = J-^s (42) 

"mm y ""mm 

where all quantities in the above depend on L\> ox . Fig- 
ure [H presents the comparison between the estimates of 
the linear bias from the simulations and the values ob- 
tained from the Shcth-Tormcn formula. The top and 
bottom panels show the results for &hm and 6hh, respec- 
tively. Again the solid red and empty blue circles denote 
the results from the cubical and spherical window func- 
tions, respectively. The Sheth-Tormen theory is repre- 
sented by the thick green dashed line. 

Considering bhm (top panel), the first thing to remark 
is that the bias appears extremely flat over all of the 
scales probed - for the mean of the realizations the bias 
relation is linear to better than 1% precision. Secondly, 
the peak-background split model of Sheth & Tormen pre- 
dicts this value astonishingly well: b = 1.498. 

Turning our attention to &hh (lower panel), the raw 
simulation measurements (upper set of points) indicate 
that on scales Lbox > 200ft.~ 1 Mpc, the bias displays 
a weak scale-dependence and is roughly ~ 3% higher 
than the Sheth-Tormen prediction. However, on smaller 
scales nonlinear effects are apparent and the overall am- 
plitude is steadily increasing with decreasing scale, be- 
ing ^ 7% higher than the Sheth-Tormen prediction for 
^box = 50/i _1 Mpc. The figure also shows the impor- 
tance of correcting a^ h for shot-noise when making es- 
timates of the bias. The upper and lower set of points 
in this panel denote the uncorrected and corrected es- 
timates, respectively. The shot-noise correction reduces 
the discrepancy between the simulations and linear the- 
ory to within ~ 2% for Lbox > 200 ft. _1 Mpc, however the 
nonlinearity on smaller scales remains. 

Both cubical and spherical window functions yield very 
similar results. In the rest of this work we shall employ 
the Sheth-Tormen bias, since on the scales of interest we 
have shown that it is at worst < 5% compared to the 
average bias of the haloes in our simulations. 

Finally, we mention that for the analytical results in 
the next sections, we shall use: (i) the volume variance 
measured from the matter-matter power spectrum with a 
cubical window function, and not the analytical variance, 
given the discrepancy seen in Figure [T] The cubical win- 
dow function is a natural choice, since our simulations 
also have this geometry; (ii) the Sheth-Tormen bias; (iii) 
the Sheth-Tormen mass function. 

5.3. An estimator for the mass function covariance 

We estimate the mass function covariance matrix from 
the ensemble of 40 simulations of the fiducial cosmolog- 
ical model, described in fj4] As we will show shortly, 
this number of realizations is insufficient for a reliable 
estimate of the covariance matrix. In order to over- 
come this problem, we have adopted the simple strategy 
of subdividing the volume associated with each realiza- 
tion into a set of smaller cubes. In particular, we di- 
vide each dimension of the original cube by 2, 3, and 4. 
Hence, each cube of 1500 3 h~ 3 Mpc 3 is partitioned into 
8, 27, and 64 subcubes with corresponding volumes of 
750 3 /i" 3 Mpc 3 , 500 3 /i- 3 Mpc 3 , and 375 3 h~ 3 Mpc 3 , re- 
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Fig. 3. — Comparison between the predicted and measured fractional error on the halo mass function as a func- 
tion of halo mass. The four panels show the results obtained when the the sample volume length is taken to be: 
Lbox = {1500 , 750 , 500 , 375 } /i _1 Mpc. In each plot, the dashed blue lines denote the fractional Poisson error; the red dot- 
dashed lines denote the pure sample variance error; and the solid lines represent the total. All errors have been rescaled to a 
total survey volume of V = 135 h~ 3 Gpc 3 . 



spectively. The 'subcubing' procedure thus provides us 
with 40, 320, 1080, and 2560 quasi-independent realiza- 
tions. We note tha t this strategy was also adopted by 
iCrocce et al.l (|2010D . who used it to compute sample- 
variance error bars on the mass function in the MICE 
simulations. However, it has never been employed to 
compute the covariance matrix of counts. 

One potential disadvantage of this approach, is that 
the realizations thus obtained are not perfectly indepen- 
dent, since there will be modes with wavelength of the 
order of the initial box size 1500 h^ 1 Mpc, which will po- 
tentially induce some covariance between the structures 
in each set of subcubes. However, as described in Ap- 
pendix [B] we have checked that this effect is of marginal 
importance. We shall therefore treat the measurements 
in each subcube as providing essentially independent in- 
formation. Conversely, the subcubing approach should 
actually be thought of as the most relevant scenario, 
since in the real Universe there is no cut-off in the power 
spectrum on scales larger than the survey. As we demon- 
strated in Figure [TJ the cut-off scale in the simulations 
dramatically affects the behaviour of the density variance 
on the scales of the box. Hence, studying the mass func- 
tion covariance using simulations that do not account for 
power on scales larger than the box modes, may in fact 
lead to incorrect inferences about the real Universe. 



Our estimator for the covariance matrix can be ex- 
pressed as follows. Let N Iuns be the total number of in- 
dependent simulations in the fiducial suite, and N sc the 
number of subcubes per simulation that we consider. For 
each subcube size, we estimate the average mass function 
as: 

N 1 Ntot 

^M a ) = — J2 > ( 43 ) 

i=l 



V sim AM a N t , 



where we defined N to t = JVruns * N sc and A^, Q is the 
number of counts in the i th subcube and mass bin a; 
V s i m = 1500 3 /i _3 Mpc 3 , and N mns = 40. We estimate 
the mass function covariance between mass bins a and 



M 



a/3 : 



n{M a )n{M p ) 



AM a AMp JV to t .4^ 



Ntot 

V Ni.„N~. 



(44) 



Note that in the above equation we subtract off the mean 
mass function averaged over all subcubes and all realiza- 
tions in bins a and (3. In order to check that the co- 
variance matrix which we present below, is not affected 
by our choice of the mean density of haloes, we recom- 
pute it using an alternative method: we determine the 
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mean density for each realization and subtract it from the 
counts in the subcubes of that realization. This alterna- 
tive is described in Appendix [Bj However, the results 
obtained from both methods are consistent. 

The covariance matrices of the counts and the mass 
function are related through the equation, 



C a0 = V*AM a AM p M a/3- 



(45) 



In the following sections we present measurements 
made at z = 0. The mass function analysis is car- 
ried out for 12 logarithmically spaced bins, going from 
(10 13 < M [H^Mq] < 10 15 ). Finally, let us make the 
clarification that when we refer to 'halo mass', we mean 
the mass returned from the FoF algorithm. 

5.4. Measurements: variance 

Figure [3] presents the fractional errors on the mass 
function, a[n(M)]/n(M), from both theory and simu- 
lations, for the subcube sizes mentioned in i)5.3l In all 
panels the points denote the measurements from the sim- 
ulations. The theoretical predictions of Eq. (f25"|) . are es- 
timated for a single realization of given size Lbox, follow- 
ing the recipe at the end of 8 35.21 Then the variance is 
rescaled by 1 /N to t , so that the fractional errors in all four 
panels correspond to a total volume of 135 h~ 3 Gpc 3 . 

The agreement between the theory and the measure- 
ments is very good, with a slight difference at the low- 
mass end for the subcubes considered. This difference 
does not occur when the estimate is made using the full 
simulation boxes to estimate the variance (see top left 
panel of Figure [3]) . We also note that the Poisson model 
(dashed lines) agrees well with the simulations at the 
high-mass end. However, at lower masses, the variance 
becomes dominated by the sample variance, as given by 
the first term of Eq. (|2"5|) . For a mass bin a the latter is 
simply: 

Wt 1 * 1 -*™' (46) 
n(M a ) 

and this is denoted in Figure [3] by the dot-dashed lines. 

On comparing all four panels, we observe that with the 
exception of the first panel with L^ ox — 1500 h~ 1 Mpc, 
the results are almost indistinguishable. This is quite in- 
teresting, since for these subcube volumes, Fig. [I] shows 
f(^i) to be a decreasing function of Lbox- For a given 
mass bin we would expect the errors for the smaller sub- 
cube measurements to be significantly larger. This is 
indeed the case, but the fact t hat w e use the variance on 
the mean, i.e. we divide by \J N to t , leads to results that 
are very similar. 

The slight difference between the measurements and 
the predictions is not easy to understand, since for the 
theoretical estimation we use cr(V^) measured from the 
simulations. This is done for all subcubes, so we do take 
into account that modes with wavelength larger than the 
subcube size may contribute to the covariance in the sub- 
cubes. The limit is set by the size of the original simula- 
tion box L S i m = 1500 /i _1 Mpc. However, the bias is the 
Sheth-Tormen prescription, which Figure [2] shows to be 
slightly lower than the one measured from the halo-halo 
power spectrum. This effect might be more pronounced 
for the small-mass bins, but more work is needed here to 



arrive to a definitive conclusion, and we defer this to a 
future study. 

Before moving on, we note that this startling agree- 
ment for the fra ctional errors on the mass function was 
noted before bv iCrocce et a l. (2010). In that work the 
variance on a given subcube scale was computed theoret- 
ically using the linear theory variance in a spherical top- 
hat taken to have the same volume as the subcube (see 
earlier discussion in Sj3]) • These authors pointed out that 
when using an ensemble of simulations with no subcub- 
ing the theory over-predicted the measurements. Here we 
have shown that there is no conflict between the theory 
and the measurements, if one uses the volume variance 
measured from the simulations. 

5.5. Measurements: covariance 

Figure|4]presents the theoretical mass function correla- 
tion matrix from Eq. (|29|) versus the measured one. The 
left panels show the predictions, obtained in the same 
way as in Figure |31 and the right ones the measurements. 
From top to bottom, the following subcube sizes are con- 
sidered: Lbox = 375, 500, 750, 1500 h^Mpc. 

The figure reveals a remarkable agreement between 
measurements and theory: the trend observed in Fig- 
ure[3]is also present here, with the predictions marginally 
larger than the measurements in some of the mass bins. 
We find that the measurements are strongly covariant: 
for clusters with M << 3 x lO 14 /i _1 M , the cross- 
correlation coefficeint is r > 0.5. Only for the highest- 
mass clusters does the covariance matrix become close 
to diagonal. The exception is for the ensemble of cubes 
with ibox = L s i m . In this case the realizations ap- 
pear to be only weakly correlated with r > 0.1, for 
M < 3 x 10 13 /i _1 M Q . However, as was discussed above 
and shown in Figure [H the behaviour of <r(V^) at the 
simulation box scale is not representative for the real 
Universe, owing to the absence of power on larger scales. 
Had we run a simulation of larger volume, then the vol- 
ume variance on the scale Lbox = 1500/i -1 Mpc would 
have been significantly larger. 

Figure[5]presents the same information as Figured) but 
in a more quantitative format. The plot has 12 panels, 
with each panel depicting a single row from the corre- 
lation matrix, i.e. rij(Mi, Mj) vs Mj, with M, fixed. 
In this plot the solid triangles denote the measured cor- 
relation coefficient, while the empty squares represent 
the theory predictions. For clarity, we show results only 
for the box sizes 1500, 750, 375 h~ 1 Mpc, represented by 
the magenta, green, and blue symbols respectively. It is 
clear from this figure too that the theory predictions and 
the measurements are in remakably good agreement. On 
comparing the correlation coefficient for different sub- 
cube sizes, we again note the similarity of these results, 
despite the variation in cr(V^): just as in Fig.|3l the co- 
variance on the mean leads to the observed similarity. 
The exception is for the L^ ox — L s - lm cubes, and we offer 
the same explantion for this as noted above. 

We conclude this section by stating that Eq. (|25|) gives 
a very reliable prediction for the mass function covari- 
ance, provided one employs the true variance within the 
volume. 

6. COSMOLOGICAL INFORMATION FROM THE MASS 
FUNCTION 
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Mass [h 'Mq] Mass [h l yi s ] 

Fig. 4. — The correlation matrix of the cluster mass function r a p, i.e. Eq. (I29[) . The left and right columns show 
the results from theory and simulations, respectively. From top to bottom, the size of the sample volume is given by: 
ibox = {375, 500, 750, 1500} ft, - Mpc. The theoretical predictions for the correlation matrix are generated using the esti- 
mate of a mm (V^) measured directly from the simulations. Note that in the bottom right panel we plot \rij\, so as to maintain 
the same heat-bar intensity scale as in the other plots. 
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L box =1500, 750, 375 h" 1 Mpc 
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Fig. 5. — Rows of the cluster mass function correlation matrix r a p as a function of the mass scale Mp, with M a fixed. Each 
of the 12 panels shows the results for one of the 12 rows of r a p. In all panels, the theoretical predictions and the measurements 
from the simulations are denoted by the empty squares and solid triangle symbols, respectively. The magenta, green, and blue 
colours represent the sample volume sizes I/box = {1500, 750, 375} /i _1 Mpc, respectively. 
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In this section we examine how the cosmological in- 
formation content of the cluster mass function changes, 
when we exchange the standard Poisson assumption for 
the more complex likelihood models of Eqs (|10[) and (fl4|) . 



6.1. Fisher information 

In all cases we shall use the standard definition of the 
Fisher information (for an excell ent review of F isher ma- 
trix techniques in cosmology see lHeavensI 2009): 

'*¥) ■ < 47 » 

OPaOPb I 

where p a and pb are elements of the cosmological model 
parameter set upon which the likelihood depends. From 
the Fisher matrix, one may obtain an estimate of the 
marginalized errors and covariances of the parameters: 



a Pa,Pb ~ ^ \paPbl 

as well as the unmarginalized errors: 

a , = rp_ . 1-1/2 



(48) 



v a - [Fp aPa ]- x ' 2 . (49) 

6.2. The Poisson Fisher matrix 

In the case of Poisson errors for each cell and mass bin, 
then using Eqs ([2]) and (fl"3j) we write: 

In£ = ^lnP(iV; >a |m Vl ) 

-m ita + N it0l lnTOj, Q - lniV JjQ !] . (50) 



On partially differentiating the above expression with re- 
spect to parameters p a and then pb, and on performing 
the ensemble average, one finds: 



F, 



'Poisson \ ^ 



PaPb 



La 



dp a dpb m i;Q 



(51) 



6.3. The Gaussian Fisher matrix 



As was shown earlier, in the case of the full likelihood 
model for the counts (c.f. Eq. ([TO])), we expect the Fisher 
matrix to be significantly modified from the Poisson case 
in the region of many counts per mass bin. In this limit, 
the likelihood is given by Eq. (TT4"f , and we have the stan- 
dard resu lt for the Fisher info rmation for a Gaussian like- 
lihood (|Tegmark et al.lll997ft : 



1. 



X"Lrauss _TY 

^PaPb ~ 2 



dp a dpb 



P S-^.(52) 

Op a Opb 



6.4. The Gauss-Poisson Fisher matrix 

LH04 developed an approximation for the Fisher ma- 
trix, which interpolates between the correct forms for 
the information in the limit of rare peaks and sample- 
variance-dominated counts. Their expression is: 
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where C = M + S, and M is a diagonal matrix with the 
elements defined in !j2l 



6.5. Estimating Fisher matrices from simulations 

In order to evaluate all of the expressions for the Fisher 
matrices presented in the previous sections, we require 
knowledge of three quantities: the partial derivatives of 
the mean counts with respect to the parameters dm./ dp a ; 
the inverse of the total covariance matrix C ; and 
the derivative of the sample variance covariance matrix 
dS/dpa- In this section we shall use numerical simula- 
tions to directly evaluate all of these quantities. 

We first measure the halo mass function for each 
of the variational cosmologies described in <j4] With 
this information we are then able to numerically ob- 
tain the derivatives dm~/dp a for the simulated param- 
eters p a € {O m , erg, n, w}. When computing the mass 
function derivatives, we reduce the effects of cosmic vari- 
ance on the estimates, using the fact that the first 4 
simulations of the fiducial cosmology have matched ini- 
tial conditions with the 4 variational-cosmologies simula- 
tions. Hence, our reduced-cosmic-variance estimator for 
the derivatives can be written: 



dn a 
dpb 



-E 



d log n a 
dpb 



(r) 



(54) 



where n a is the average mass function for mass bin a, 
estimated from all 40 independent realizations; -/V var = 4 
is the number of the variational simulations; r denotes 

_ (r) 

the simulation realization going from 1 to N V&T ; n a is 
the mass function in the fiducial case, estimated for each 
of the 4 realizations that have matched initial conditions 
to the variational-cosmol ogie s realizations (for an explict 
defintion of this see Eq. (|B1[) ). The logarithmic deriva- 
tives are estimated as: 



d log n a 
dpb 



(r) 



(Pb + A 6 ) 



(pb - A 6 ) 



2A b n a r \p b ) 



(55) 



Note that since we estimate drn/dp a using double-sided 
derivatives, we may take larger step sizes in the pa- 
rameters to compute the derivat ives than would be a l- 
lowed for single sided derivatives (|Eisenstein et al.f l999'). 
For the former case, the errors in the derivatives are 
of quadratic order in the step size: i.e. A[dm./dp a ] ~ 
(Ap a ) 2 d 3 m/dp a /6. Thus parameter step sizes of 20% 
and 10% should correspond to relative errors of roughly 
4% and 1% in the derivatives, respectively. In actuallity, 
the true accuracy of the derivatives also depends on the 
value of the third partial derivative. 

In Figure [S] we show simulation measurements of the 
average mass functions for the fiducial and variational 
cosmologies. This figure makes very clear not only the 
sensitivity of the mass function to the cosmological pa- 
rameters considered, but also the halo mass range over 
which most of it occurs. Changes in fi m and the slope of 
the primordial power spectrum n impact the mass func- 
tion for the whole range of halo masses. The low-mass 
end is less sensitive to variations in as, while the dark 
energy equation of state parameter w barely affects the 
mass function. 

In the smaller panels of Figure |5] we show the deriva- 
tives of the halo abundance, estimated using Eq. (|55]h 
The error bars are computed as errors on the mean of 
the iV var = 4 realizations, as they are also for the mass 
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Mass [h _1 M ] 



Fig. 6. — Top section of each panel: Dependence of the z = cluster m ass function on cosmolo gy, as a function of cluster 
mass. Symbols denote measurements from the simulations and lines depict the Sh eth fc Torm en (1999) mass function. The green 
colour represents the fiducial model, whereas the red/blue colours are for the plus/minus variations in the parameters. Bottom 
sections: Logarithmic derivatives of the cluster number counts with respect to the considered parameters, as a function of 
cluster mass. Points with errors denote measurements from the simulations (c.f. Eg. 1551) . the error bars being on the mean. 
Lines denote the Sheth- Tormen predictions. 



functions in the larger panels. The f2 m -derivative is al- 
most constant and large for all bins, while the one 
monotonically increases from at the low mass end to 
a large value at the high mass end. The behaviour of 
the spectral-index-derivative is quite interesting, as it 
changes sign at M ~ 3 x 10 14 /i _1 M Q and becomes neg- 
ative at the high mass end. Its overall variation is not as 
large as in the case of O m and as, which will be better 
constrained by the halo abundance. 

Another interesting finding of this exploration concerns 
the ui-derivative, which should be at redshift accord- 
ing to linear theory and the Sheth- Tormen mass func- 
tion. We find it to be small and positive, ~ 0.05, for 
most of the mass range considered, and rising slightly 
to ~ 0.1 at the low-mass end. The w-derivative does 
not vanish because in reality the mass function depends 
on the full nonlinear growth history. This encompasses 



the growth of structure at all redshifts, and thus makes 
the present day halo abundance sensitive to w. These 
results are consistent with the findings in ear lier studies 
(jLinder fc Jenkindl2003t Pennings et al.ll2010f ). 

We next follow the recipe of §5.3| to estimate the co- 
variance matrices in each of the variational cosmological 
models. From these estimates we are then able to form 
the partial derivatives of the covariance with respect to 
the cosmological parameters: dC/dp a . Again, as was 
done for dm./dp a , we take advantage of the matched in- 
tial conditions to reduce the cosmic variance on the esti- 
mates of the partial derivatives of the covariance matrix. 

6.6. Forecasted errors 

Having obtained all of the necessary ingredients we are 
now in a position to evaluate the Fisher information di- 
rectly from the simulations. 
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Fig. 7. — Fractional Fisher-matrix errors on the cosmological parameters, as a function of the minimum cluster mass used. 
The four panels show the results for the cosmological parameters: p a £ {fi m , erg, n, w}. In all panels, the symbols show the 
estimates made from the iV-body simulations, with varying assumptions about the form of the cluster likelihood function. Solid 
green squares denote the Poisson errors obtained with Eq. (I51|) ; solid red circles denote the errors obtained from the second term 
of Eq. (|53[) ; blue triangular shaped symbols denote the errors derived from the trace-term in Eq. (I53|) . 



Figure [7] shows the cumulative fractional Fisher errors, 
Ap a /pa d , estimated using Eq. (j49|) . as a function of the 
minimum cluster mass, and for the four cosmological pa- 
rameters that we consider. The results obtained for the 
various subcube sizes are almost identical with the excep- 
tion of the case where Lbox = ^sim^ as explained earlier, 
the underestimate of the variance on scales of the sim- 
ulation makes the estimate of the mass function covari- 
ance, and hence the Fisher errors unrealistic. For brevity 
we shall present only the findings for L = 375 h~ 1 Mpc, 
which we consider very reliable. 

For our fiducial survey, we adopt parameters rele- 
vant for fu ture all-sky X-ray cluster surveys, such as 
eROSITA (|Predehl et all 120101 ) . This mission will be 
able to target intermediate mass range clusters, and not 
just the most massive objects in the Universe as is the 
case for current and past surveys. We adopt a total 
survey volume of V ~ 13.5 h~ 3 Gpc 3 , and we rescale 
our measured covariance matrices to this volume. For 
this comoving volume at z = 0, we find in the simu- 
lations approximately 4.5 x 10 6 halos in the mass in- 
terval [1,5] x 10 13 /i _1 Afo, 5.4 x 10 5 halos in the inter- 
val [0.5, 1] x 10 u h- 1 M @ , 2.3 x 10 5 halos in the interval 
[1, 6.5] x 10 14 h~ 1 M Q , and 8000 halos with masses larger 
than the latter limit. 

In Fig. [7] the solid green squares denote the results ob- 
tained for the Poisson Fisher matrix, as given by Eq. (fSTj) . 



The solid red circles denote the errors resulting from 
the second term of Eq. ([53| . The blue triangular-shaped 
symbols denote the errors obtained from only the trace- 
part of Eq. (|53"|) . where instead of S, we have used the 
covariance matrix from our simulations C (c.f. Eq. (|45p). 
We do not expect that replacing C with S will change 
our conclusions concerning the information carried by 
this term, except to possibly make the errors larger. 

As expected, for all of the cosmological parameters 
considered, the fractional errors obtained from the Pois- 
son approximation are smallest. Including the full covari- 
ance matrix, as in the second term of Eq. (|53[) , reduces 
the amount of information, and this results in a signif- 
icant increase in the fractional errors. For the case of 
M m i n ~ 10 13 /i _1 M Q , the errors are roughly a factor of 
^3 larger when the full-covariance is used as opposed to 
the Poisson case. When M m - m ~ 10 14 Ii~ 1 Mq, the errors 
are only a factor of ^2 worse. For the rarest objects, 
where the covariance becomes almost diagonal, the er- 
rors from the two methods are very similar. We find 
that the trace part of Eq. (1531 contributes negligibly to 
the information, and if this term is taken separately, it 
yields errors that are roughly one order of magnitude 
larger than those from the second term. 

Let us explore the consequences of this last result a lit- 
tle further. Consider the Fisher matrix given by Eq. ([53"1) . 
if the first term on the right-hand-side is negligable, then 
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the information about each cosmological parameter en- 
ters the system only through the derivatives of the model 
mean with respect to the parameters. Since the model 
here is the mean counts, the bias provides no information. 
However, the amplitude of the elements of the informa- 
tion matrix can be modulted by the inverse covariance 
matrix. Owing to the fact that increasing the elements 
of the covariance matrix only leads to a smaller inverse 
covariance, we thus conclude that, adding the variance 
from the bias can only ever decr ease the Fisher inf orma- 
tion. However, as discussed in iLima fe Hul ([20051 ). the 
information content of the first term of Eq. (|53|) , becomes 
of great importance in the presence of a scatter between 
the true and observed mass. 

Note also that the cumulative dependence of the errors 
on the mass bins, can partly be understood by examining 
the behaviour of the derivatives as shown in Fig.[6l The 
errors flatten out at those points in the mass range where 
the derivatives of the parameters are close to 0, as in 
the case of as at the low-mass end, or n at masses ~ 
3 x IO^/T^Mq. 

Finally, we emphasize that the forecasts that we make 
above are to illustrate the importance of going beyond 
the Possion likelihood approximation and should not be 
taken as serious predictions for a potential survey. The 
cosmological dependence that we have considered here 
arises strictly from the mass function. In order to make 
a realistic forecast we would have to take into account a 
number of observational factors: realistic survey geome- 
tries; the evolution of the mass function with redshift; 
the evolution of the volume element with the cosmologi- 
cal model; and the evolution in the minimum detectable 
mass at each redshift; and a scatter in the relation be- 
twee n the observed mass proxy and the true cluster mass 
(see iMarian fe Bernstein! 120061 for an example of fore- 
casting weak lensing cluster counts.). 

7. SUMMARY AND CONCLUSIONS 

In this paper, we have studied the covariance of the 
halo mass function, and the cosmological information 
content of such data. We adopted a two-line attack on 
these problems: the first line was theoretical and we de- 
veloped an analytic model to explore these issues; the 
second was the use of a large ensemble of numerical sim- 
ulations to measure directly all quantities of interest. 

In j|2l we summarized the counts - in-cells formalism 
([Peebles! fl980l IHu fc Kravtsovl 120031 ). and developed it 
for application to deal with cluster counts in multi- 
ple mass bins. We described the Gauss-Poisson likeli- 
hood function for the counts in cells with multiple mass 
bins. The express ion was analogous to that derived by 
ILima fc Hul (|2004[ ) for multiple cells and a single mass 
bin. 

In we used this framework to derive a formal ex- 
pression for the covariance of the halo mass function and 
the cross-correlation coefficient, i.e. Eqs ([231) and (|2"§)) , 
respectively. We found that there were two terms con- 
tributing: a Poisson shot-noise term, which dominates 
in the limit of rare clusters; and a term associated with 
the sample variance, which is dominant for abundant 
clusters. This express ion is analogous to the results of 
IHu fc Kravtsovi (|2003l) for multiple cells and a single 
mass bin. The expression was found to depend on three 
quantities: the cluster mass function; the cluster bias; 



and the variance in the survey volume. 

In SjH we presented the details of our large ensemble of 
numerical simulations: 40 simulations of a fiducial model 
and 32 simulations of modified cosmological models. 

In <j5j we made a rigorous comparison of the results 
from the theoretical modelling with those obtained di- 
rectly from the numerical simulations. We measured the 
variance of matter and cluster fluctuations in cells of var- 
ious sizes and found, for spherical and cubical top-hat 
cells, that the simulations and theory predictions were 
discrepant for large cell sizes. We showed that this was 
entirely attributable to the difference between the dis- 
crete lattice structure of the Fourier space in the simula- 
tions, and the continuum of Fourier modes in the theory 
integrals. The cubical and spherical top-hat simulation 
results were in good agreement, except on the largest 
scales where simulation box-scale effects were important. 

We also measured the halo bias in cells of various sizes 
from the simulations. We found that the bias from the 
halo-mass cross-variance showed very little scale depen- 
dence over the range L hox = [50, 1500] /i _1 Mpc, whereas 
that from the halo auto- variance showed significant scale 
dependence, before and after the shot-noise correction. 
We found that the iSheth fc Tormenl (|1999f ) model was 
an excellent fit to the former and a reasonable fit to the 
latter. 

We then measured the covariance of the mass function 
in the simulations. To increase the number of realiza- 
tions, we used the strategy of subdividing each large sim- 
ulation volume into a set of smaller subcubes. We found 
that the estimated covariances were in excellent agree- 
ment with the theoretical predictions. This was under 
the condition that we used the actual variance of mass 
fluctuations measured in the simulations. 

In Sj6j we employed the Fisher matrix formalism to ex- 
plore the information content of the cluster counts. Us- 
ing the more realistic likelihood functions, we demon- 
strated numerically that the Poisson likelihood model 
only provides a reasonably accurate description of the 
data for clusters that are more massive than M > 5.0 x 
1O 14 /i _1 M0. Future surveys that aim to target cluster 
samples with masses M < 5 x 10 14 /i _1 M , must adopt 
mo re sophistica t ed lik el ihood analys i s, suc h as discussed 
by ILima fc Hul (|200l . IHu fc Cohnl (1200(1 and here in, 
which take into account the full covariance matrix of the 
counts. Otherwise, significant underestimates of the true 
errors will occur. 

There are a number of possible future directions for the 
wo rk that we have pre sented here. Firstly, as pointed out 
by ILima fc Hul ([2005D , one of the main uses of adopt- 
ing the counts-in-cells approach is that it helps to lift 
the degeneracy between nuisance parameters, which are 
involved in calibrating the cluster masses, and the cos- 
mological parameters. This occurs beacuse the sample 
varaince depends on the bias of the clusters, which has 
a different behaviour with cosmological parameters than 
the counts. Whilst we have shown explicitly that the 
terms in the Fisher matrix that depend on the derivatives 
of the covariance matrix, and hence derivatives of the 
bias, do not carry a great deal of cosmological informa- 
tion by themselves, it will be interesting to see whether 
for a more realistic scenario, where one must marginalize 
over these nuisance parameters, the self-calibration can 
be successfully performed to restore the lost information. 
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We also note that the counts-in-cells technique has 
been highlighted as a power ful means fo r constraining 
prim o rdial non-Gaussian ity (|Qguril 120091 iCunha et al.l 
[20101 iMarian et afll20Tl . It is of some importance to 
explore this approach using numerical simulations, since 
it is not clear whether the extension of the current for- 
malism to such modified cosmological models works in 
practice. 
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Fig. 8. — The r.m.s. density V3,ricLllCG (7 mm (V^i) as a function of the sample volume size Lbox- In each panel, blue empty and 
solid red circles denote theoretical predictions made on the Fourier space lattice using the spherical and cubical top-hat filter 
functions, respectively. The predictions made using a continuum of Fourier space modes are denoted by the dashed blue and 
solid red lines, respectively. The top panel compares the results when the simulation box size is taken to be Lsim = 1500 h~ 1 Mpc. 
The bottome panel shows the same but for the case where L s i m = 6000 /i _1 Mpc. 

APPENDIX 

THE VOLUME VARIANCE 

In this section we investigate the impact of systematic effects on the volume variance, which arise due to the 
anisotropic lattice in Fourier space and also low- and high-A; truncation of the matter power spectrum. 

Fourier lattice versus a continuum of modes 

As was described in i j3.4l the matter variance in the volume is a key quantity for correctly evaluating the covariance 
matrix of the cluster counts. Also, as was shown in H5.ll there is a discrepancy between the theoretical predictions 
and meaurements in simulations obtained for er(V^). We now investigate the origin of these discrepancys. 

We start by examining the importance of the discrete cubical Fourier lattice, which is used in the estimates from 
the simulations, and the continuum of Fourier modes, which is used to evaluate the theory. We start by generating 
a Fourier lattice as in the simulations, where each lattice point is spaced from the next one, along each dimension, 
by fcf = 27r/L s ; m . Then at each lattice point, we compute the magnitude of the k- vector and evaluate the linear 
theory power spectrum at that point. er(VJj) is then obtained as described in Eq. ()39[) . by summing up the grid of 
power spectra values multiplied by the square of the appropriate window function. The top panel of Figure [5] shows 
the results of this exercise for both the spherical and cubical window functions. We also compare this to the results 
obtained from the theory, assuming a continuum of modes. The results that we find for the theory evaluated on the 
cubical mesh, are in remarkably good agreement with the measurements from the simulations presented in Fig.[TJ 

To be sure that the discrepancy is due to the lattice, we should expect that as the simulation box size becomes 
significantly larger, the results for the lattice should approach those of the continuum. We test this by regenerating 
the Fourier lattice, but this time taking L s ; m = 6000 /i _1 Mpc, and keeping the maximum Fourier mode the same as 
before. The results of this exercise are shown in the bottom panel of Fig. [8] We clearly see that the results are now in 
much better agreement and for the same cell sizes as in the upper panel. 

Thus we are led to conclude that in matching the results from the simulations we must be mindful to take into 
account the anisotropic lattice structure of the Fourier space to obtain accurate comparisons between the theory and 
the simulations. This then further justifies our choice of using the cr(V^) measured in the simulations to make the 
predictions for the covariance of the counts. 

Finally, these results also act as a cauationary tale: when interpreting the results of numerical simulations on very 
large scales, we should take more care in asigning the power to the lattice cells in the initial conditions. We should use 
methods that supress this discretization. For instance, it would seem more sensible to compute the power averaged 
over a lattice cell and not simply the po wer at the l attice cell point. Also including the missing zero modes may be a 
more realistic stratergy (see for exa mple Sirko 2005). 

Cut-off scales 

In Figured] we evaluated the integrals in Eqs (|35]l and ([37]) . keeping the upper and lower bounds fixed at the values 
fcmin = 2"7r/L s i m = 0.004 h Mpc -1 and fc max = 27r/i so ft = 100 /iMpc -1 . This was done for a fair comparison with our 
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Fig. 9. — Dependence of the r.m.s. density variance on the lower limit fc m in of the fc-space integrals, as a function of the 
sample volume size Lbox- The thin solid red, dashed green and dot-dashed magenta lines denote results obtained for a cubical 
top-hat window function, where fc m i n = 27r/L s i m with Z/ S i m = {750, 1500, 3000} /i -1 Mpc, respectively. The thick solid, dashed 
and dot-dashed lines represent the same, but for the case where the filter function is a spherical top-hat. The thick and thin 
blue dotted lines denote the same as above except this time the lower limit of the fc-space integrals is given by fc m i n = 2ir/Lbox- 



simulations, which do not have modes larger than the simulation box L s i m = 1500 /i _1 Mpc, nor structures smaller 
than the softening scale, Z so f t = 0.06 ft, Mpc. We now present a short discussion of how the mass-fluctuations- variance 
<7(Vfj_) depends on the cell volume and the cut-off scales in the power spectrum. 

For the large cell sizes that we arc interested in, i.e. Lbox > 50/i Mpc, we find no dependence of cr(V^) on fc max , 
for the range of values fc max = [1, 100] /iMpc -1 . 

For the lower cut-off scale fc min , the situation appears to be more complex. In Figure |H1 we show the result of 
computing the mass-fluctuations-variance averaged in cubical and spherical top-hat volumes, as a function of the 
cubical cell volume (recall that we take the radius of the spherical top-hat cell to be R = (3/47r) 1//3 Lbox)- In the plot 
we consider the values of er(V^) for three different simulation sizes: L sim = {750, 1500, 3000 /i _1 Mpc}. These box sizes 
correspond to the: solid red, long-dashed green, and dot-dashed magenta lines, respectively. The thicker/thinner lines 
in the plot depict the spherical/cubical top-hat smoothing. 

For the case of the spherical top-hat filter, we find that the variance for fc m j n = 27r/750 h Mpc - = 0.008 h Mpc~ , is 
roughly a factor of ^2 times smaller than the variance obtained when fc m i„ = 27r/3000 = 0.002 hMpc^ 1 . However, for 
the case of the cubical top-hat window function, we find that the difference in the variance for these same two values 
of fcmin, is more than an order of magnitude. 

In Figure^ the thick dotted blue curve presents predictions for cr(V^) with the spherical window function, but where 
we now take the lower limit fc m ; n = 2ir/L. The thick and thin dashed blue lines show the same, but for the case of 
the cubical filter function. For this case, the thin line is obtained when the linear theory matter power spectrum is 
used, and the thicker line shows the results obtained when the nonlinear power spectrum from halof it flSmith et al.1 
120031 ) is employed. The differences are very small. Thus using the linear theory power spectrum for the mass variance 
is quite reasonable on these scales. The main point of this last example, is to show that for large cell sizes, the results 
for cr(V^) are very sensitive to the presence/absence of power on very large scales. 

CONVERGENCE OF THE COVARIANCE MATRIX 

Covariances from individual simulations 

Here we consider an alternate approach to estimating the covariance of the cluster counts. We are concerned that, 
if there is a significant variance of the cluster counts on the scales of the simulation cube, then by computing the 
covariance around the mean cluster mass function averaged over all simulations, we are overestimating the covariance. 
To anwer this question, we adopt the stratergy of using the sub-cubes in a single simulation to make an estimate of 
the covariance, and finally we then average these estimates over all the simulations. 
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Fig. 10. — The chessboard test: we compare the mass function correlation matrix measured from 'white' and 'black' subcubes- 
see the Appendix text. The results are very similar for both subcube sizes considered, 375 3 ft -3 Mpc 3 and 250 3 /i _3 Mpc 3 . 
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For each simulation run we therefore have: 
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where is the number of counts in the r th run, i th subcube and a th mass bin. The covariance for each run is: 
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We have checked that using Eqs (1B1[) and (|B2I) does not change the measured mass function covariance in any significant 
way. We therefore conclude that the method of estimating the covariance described in is not biased by the 

estimates of the mean density. 

The chessboard test 

When dividing a big simulation box into smaller subcubes, the largest wavelength modes may affect the observables 
measured in the smaller cubes. In the case of clusters, some of the subcubes may have very different mean counts than 
others, and in general, the smaller the sub-boxes, the larger the expected covariance between them. This is also true 
for real surveys, which measure observables in a finite volume of the Universe: some of these observables are impacted 
by modes larger than the size of the survey. 

In order to check the validity of our approach, we measure the covariance of the mass function using subcubes that 
are not adjacent, and should therefore be less covariant. We shall refer to this as the 'chessboard test', as its 2D 
analogue would be similar to using only the white or the black squares of a chessboard to compute the mass function 
covariance. This test has the limitation that large mode correlations can span more than just 2 subcubes, particularly 
if the latter are small. Nevertheless, if we find the covariance measured from the 'white' subcubes different from that 
obtained from the 'black' ones and also different from the all-subcubes-covariance, then our box-division method is 
flawed. 

We perform this test for the conservative values n = 4, 6, i.e. we consider 4 3 and 6 3 subcubes, with volumes 
375 3 h~ 3 Mpc 3 and 250 3 /i~ 3 Mpc 3 , respectively. The result is shown in Figure [TU] There is no major difference 
between the 'white', 'black', and total mass function correlation matrix (c.f. Fig. [4]). We conclude that modes with 
wavelengths smaller than the size of the subcubes considered here, are not explicitly responsible for generating the 
mass function covariance. 



